#Create Vectors for coefs and standard errors for each model, and variable names
    #note that we  exclude "margin squared since it doesn't appear in either model

coef.vec.1<- c(0.08, 0.18, 0.33)
se.vec.1 <-  c(0.05, 0.04, 0.04)
coef.vec.2 <-  c(0.78, 1.74, 2.81)
se.vec.2 <- c(0.04, 0.04, 0.04)
var.names <- c("Not Very Satisfied", "Somewhat Satisfied", "Very Satisfied")

y.axis <- length(var.names):1#create indicator for y.axis, descending so that R orders vars from top to bottom on y-axis
adjust <- .025 #create object that we will use to adjust points and lines up and down to distinguish between models
    
#pdf("pekkanen_fig.pdf", height = 8, width = 7) #open pdf device
#png("pekkanen_fig.png", height = 700, width = 600)


layout(matrix(c(2,1),1,2), #in order to add variable categories and braces to left side of plot, 
    widths = c(0.1, 5))#we use layout command, create a small second panel on left side.
                #using c(2,1) in matrix command tells R to create right panel 1st
#layout.show(2) #can use this command to check results of layout command (but it must be commented out when creating PDF).

par(mar=c(8,8,.5,.5), lheight = .8)#set margins for regression plot
plot(coef.vec.1, y.axis+adjust, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.2, #plot model 1 coefs using black points (pch = 19, default = black), adding the "adjust amount" to the y.axis indicator to move points up
    xlim = c(min((coef.vec.1-qnorm(.975)*se.vec.1 -.2), (coef.vec.2-qnorm(.975)*se.vec.2 -.2), na.rm = T), #set xlims at mins and maximums (from both models) of confidence intervals, plus .1 to leave room at ends of plots
    max((coef.vec.1+qnorm(.975)*se.vec.1 +.2), (coef.vec.2+qnorm(.975)*se.vec.2 +.2), na.rm = T)),  #use na.rm=T since vectors have missing values
    ylim = c(min(y.axis), max(y.axis)), main = "")

axis(1,at = seq(-.25,3, by = .25), label = seq(-.25,3, by = .25), pretty, mgp = c(3,1,0), cex.axis = 0.8)
title(xlab="Approval")
#add x-axis and labels; "pretty" creates a sequence of  equally spaced nice values that cover the range of the values in 'x'-- in this case, integers
axis(2, at = y.axis, label = var.names, las = 1, tick = T, cex.axis =0.8)#add y-axis and labels; las = 1 makes labels perpendicular to y-axis
#axis(3,pretty(coef.vec.1, 3))#same as x-axis, but on top axis 
abline(h = y.axis, lty = 2, lwd = .5, col = "light grey")#draw light dotted line at each variable for dotplot effect
#box(bty="l")#draw box around plot
segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis+adjust, coef.vec.1+qnorm(.975)*se.vec.1, y.axis+adjust, lwd =  1.3)#draw lines connecting 95% confidence intervals

abline(v=0, lty = 2) # draw dotted line through 0 for reference line for null significance hypothesis testing

#add 2nd model
    #because we are using white points and do want the lines to go "through" points rather than over them
        #we draw lines first and the overlay points
segments(coef.vec.2-qnorm(.975)*se.vec.2, y.axis-adjust, coef.vec.2+qnorm(.975)*se.vec.2, y.axis-adjust, lwd =  1.3)#draw lines connecting 95% confidence intervals

points(coef.vec.2, y.axis-adjust, pch = 21, cex = 1.2, bg = "white" ) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color

#add legend (manually) to identify which dots denote model 1 and which denote model 2
legend("topright", c("Approval of Congress", "Approval of MC"), pch = c(19,21), cex=0.8)


#dev.off()